"""
PETSc Backend Integration

This module provides integration with PETSc.jl for high-performance
parallel linear algebra operations. Following the pattern from
LidDrivenCavity.jl.
"""

if HAS_PETSC
    using PETSc_jll
    import PETSc
    using Printf
    
    """
        PETScNSSolver <: NSLinearSolver
        
    PETSc-based linear solver for Navier-Stokes equations.
    """
    mutable struct PETScNSSolver <: NSLinearSolver
        matrix::PETSc.Mat
        ksp::PETSc.KSP
        pc::PETSc.PC
        solver_type::Symbol
        max_iter::Int
        tolerance::Float64
        
        function PETScNSSolver(matrix::AbstractMatrix, solver_type::Symbol=:gmres;
                              max_iter::Int=1000, tolerance::Float64=1e-6)
            # Convert Julia matrix to PETSc Mat
            petsc_mat = convert_to_petsc_matrix(matrix)
            
            # Create KSP solver
            ksp = PETSc.KSP()
            PETSc.KSPCreate(PETSc.PETSC_COMM_WORLD, ksp)
            PETSc.KSPSetOperators(ksp, petsc_mat, petsc_mat)
            
            # Set solver type
            if solver_type == :gmres
                PETSc.KSPSetType(ksp, PETSc.KSPGMRES)
            elseif solver_type == :cg
                PETSc.KSPSetType(ksp, PETSc.KSPCG)
            elseif solver_type == :bicgstab
                PETSc.KSPSetType(ksp, PETSc.KSPBCGS)
            end
            
            # Get preconditioner
            pc = PETSc.PC()
            PETSc.KSPGetPC(ksp, pc)
            PETSc.PCSetType(pc, PETSc.PCILU)
            
            # Set tolerances
            PETSc.KSPSetTolerances(ksp, tolerance, PETSc.PETSC_DEFAULT, 
                                  PETSc.PETSC_DEFAULT, max_iter)
            
            PETSc.KSPSetFromOptions(ksp)
            PETSc.KSPSetUp(ksp)
            
            new(petsc_mat, ksp, pc, solver_type, max_iter, tolerance)
        end
    end
    
    """
        convert_to_petsc_matrix(matrix::AbstractMatrix) -> PETSc.Mat
        
    Convert Julia matrix to PETSc format.
    """
    function convert_to_petsc_matrix(matrix::AbstractMatrix)
        m, n = size(matrix)
        
        # Create PETSc matrix
        petsc_mat = PETSc.Mat()
        PETSc.MatCreate(PETSc.PETSC_COMM_WORLD, petsc_mat)
        PETSc.MatSetSizes(petsc_mat, m, n, PETSc.PETSC_DECIDE, PETSc.PETSC_DECIDE)
        PETSc.MatSetFromOptions(petsc_mat)
        PETSc.MatSetUp(petsc_mat)
        
        # Fill matrix values
        if issparse(matrix)
            rows, cols, vals = findnz(matrix)
            for (i, j, v) in zip(rows.-1, cols.-1, vals)  # PETSc uses 0-based indexing
                PETSc.MatSetValue(petsc_mat, i, j, v, PETSc.INSERT_VALUES)
            end
        else
            for i in 1:m, j in 1:n
                if matrix[i,j] != 0
                    PETSc.MatSetValue(petsc_mat, i-1, j-1, matrix[i,j], PETSc.INSERT_VALUES)
                end
            end
        end
        
        PETSc.MatAssemblyBegin(petsc_mat, PETSc.MAT_FINAL_ASSEMBLY)
        PETSc.MatAssemblyEnd(petsc_mat, PETSc.MAT_FINAL_ASSEMBLY)
        
        return petsc_mat
    end
    
    """
        solve_linear_system!(solver::PETScNSSolver, rhs::AbstractArray, 
                            initial_guess::AbstractArray, equation_type::Symbol, 
                            verbose::Bool=false) -> AbstractArray
                            
    Solve linear system using PETSc solver.
    """
    function solve_linear_system!(solver::PETScNSSolver, rhs::AbstractArray, 
                                  initial_guess::AbstractArray, equation_type::Symbol, 
                                  verbose::Bool=false)
        
        # Convert to PETSc vectors
        b = convert_to_petsc_vector(vec(rhs))
        x = convert_to_petsc_vector(vec(initial_guess))
        
        if verbose
            println("  Solving $(equation_type) system using PETSc: $(length(vec(rhs))) DOFs")
        end
        
        start_time = time()
        
        # Solve system
        PETSc.KSPSolve(solver.ksp, b, x)
        
        solve_time = time() - start_time
        
        # Get convergence information
        reason = PETSc.KSPGetConvergedReason(solver.ksp)
        iterations = PETSc.KSPGetIterationNumber(solver.ksp)
        
        if verbose
            residual_norm = PETSc.KSPGetResidualNorm(solver.ksp)
            println("  PETSc solver converged: $(reason >= 0)")
            println("  Iterations: $iterations")
            println("  Residual: $(residual_norm)")
            println("  Solve time: $(solve_time) seconds")
        end
        
        # Convert back to Julia array
        solution_vec = convert_from_petsc_vector(x)
        solution = reshape(solution_vec, size(initial_guess))
        
        # Cleanup
        PETSc.VecDestroy(b)
        PETSc.VecDestroy(x)
        
        return solution
    end
    
    """
        convert_to_petsc_vector(vec::Vector) -> PETSc.Vec
        
    Convert Julia vector to PETSc format.
    """
    function convert_to_petsc_vector(vec::Vector)
        n = length(vec)
        petsc_vec = PETSc.Vec()
        PETSc.VecCreate(PETSc.PETSC_COMM_WORLD, petsc_vec)
        PETSc.VecSetSizes(petsc_vec, n, PETSc.PETSC_DECIDE)
        PETSc.VecSetFromOptions(petsc_vec)
        
        for (i, v) in enumerate(vec)
            PETSc.VecSetValue(petsc_vec, i-1, v, PETSc.INSERT_VALUES)
        end
        
        PETSc.VecAssemblyBegin(petsc_vec)
        PETSc.VecAssemblyEnd(petsc_vec)
        
        return petsc_vec
    end
    
    """
        convert_from_petsc_vector(petsc_vec::PETSc.Vec) -> Vector
        
    Convert PETSc vector to Julia format.
    """
    function convert_from_petsc_vector(petsc_vec::PETSc.Vec)
        n = PETSc.VecGetSize(petsc_vec)
        result = zeros(Float64, n)
        
        for i in 0:n-1
            result[i+1] = PETSc.VecGetValue(petsc_vec, i)
        end
        
        return result
    end
    
    """
        create_petsc_ns_solver(matrix::AbstractMatrix, equation_type::Symbol, verbose::Bool)
        
    Create PETSc linear solver for Navier-Stokes equations.
    """
    function create_petsc_ns_solver(matrix::AbstractMatrix, equation_type::Symbol, verbose::Bool)
        if equation_type == :velocity
            solver = PETScNSSolver(matrix, :gmres; max_iter=1000, tolerance=1e-8)
        else  # pressure
            solver = PETScNSSolver(matrix, :gmres; max_iter=2000, tolerance=1e-10)
        end
        
        if verbose
            println("  Created PETSc $(solver.solver_type) solver")
            println("  Max iterations: $(solver.max_iter)")
            println("  Tolerance: $(solver.tolerance)")
        end
        
        return solver
    end
    
    @info "PETSc backend integration loaded successfully"
else
    @info "PETSc backend not available, using Julia fallbacks"
end